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We investigate the time average mean square displacement 5^{x{t)) = + A) — 

x{t'y^dt' /{t — A) for fractional Brownian and Langevin motion. Unlike the previously investi- 
gated continuous time random walk model 5^ converges to the ensemble average {x^) ~ t^^ in the 
long measurement time limit. The convergence to ergodic behavior is however slow, and surpris- 
ingly the Hurst exponent iif = 3/4 marks the critical point of the speed of convergence. When 
H < 3/4, the ergodicity breaking parameter EB = Ysir(P)/(Pf ~ k{H) ■ A ■ when H = 3/4, 
EB ~ (9/16)(lnt) • A • and when 3/4 < H < 1,EB ~ k{H)A'^-'^"t'^"~'^. In the ballistic limit 
7? — > 1 ergodicity is broken and EB ~ 2. The critical point iif = 3/4 is marked by the divergence of 
the coefficient k{H). Fractional Brownian motion as a model for recent experiments of sub-diffusion 
of mRNA in the cell is briefly discussed and comparison with the continuous time random walk 
model is made. 

PACS numbers: 02.50.-r, 05.30.Pr, 05.40.-a, 05.10.Gg 
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I. INTRODUCTION 

Fractional calculus, e.g. d^/^/di^/^, is a powerful 
mathematical tool for the investigation of physical and 
biological phenomena with long-range correlations or 
long memory [ij. For example fractional calculus de- 
scribes the mechanical memory of viscoelastic materials 
An important application of fractional calculus is in 
the stochastic modeling of anomalous diffusion. Frac- 
tional Fokker-Planck equations describe the long time 
behavior of the continuous time random walk (CTRW) 
model, when waiting times and/or jump lengths have 
power-law distributions [E Si, 111- A different stochas- 
tic approach to anomalous diffusion is based on fractional 
Brownian motion (fBM) [6| which is related to recently 
investigated fractional Langevin equations (see details 
below) 0,ii. 

Recent single particle tracking of mRNA molecules fToj 
and lipid granules [11] in living cells revealed that time 
averaged mean square displacement (defined below 
more precisely) of individual particles remains a random 
variable while indicating that the particle motion is sub- 
diffusive. This means that the time averages are not iden- 
tical to ensemble averages. Such breaking of ergodicity 
was investigated within the sub-diffusive CTRW model 
[13, [13] ■ It was shown that transport and diffusion con- 
stants extracted from single particle trajectories remain 
random variables, even in the long measurement time 
limit. For a non-technical point of view on this prob- 
lem see [l3]- Here we consider three stochastic models 
of anomalous diffusion: fBM and the under-damped and 
the over-damped fractional Langevin equation. Except 
for limiting cases (i.e. ballistic diffusion) we find that the 
time average 6'^, in the long measurement time limit, is 
identical to the ensemble average (x^), indicating that 
these models are ergodic. Note however, that experi- 
ments on anomalous dynamics of particles in the cell are 
always conducted for finite times (due to the life time of 



the cell). Here we find the finite time corrections to er- 
godic behavior, namely we give estimates on how far will 
a finite time measurement of anomalous diffusion devi- 
ate from the ensemble average. Since the convergence 
to ergodic behavior is slow our results seem particularly 
important to finite time experiments. The problem of es- 
timating diffusion constants from single particle tracking, 
for normal diffusion, is already well investigated [l3] • 

In recent years there was much interest in non- 
ergodicity of anomalous diffusion processes. A well inves- 
tigated system are blinking quantum dots [l^ . [iTj , which 
exhibit a Levy walk type of dynamics (a super-diffusive 
process). Very general formula for the distribution of 
time ave rag es for weakly non ergodic systems was de- 
rived in jl8|, and this framework was shown to describe 
the sub-diffusive CTRW [13] . Bao et al have investigated 
ergodicity breaking for stochastic dynamics described by 
the generalized Langevin equation [20| . They considered 
the time averaged velocity variance. The latter converges 
to ksT/m in thermal equilibrium if the process is er- 
godic. It was shown [20| that under certain conditions 
the generalized Langevin equation is non ergodic (see 
also [21, 22, 23, 24, 25j). Our work, following the recent 
experiments [lO|, [ll|, considers the time average of the 
mean square displacements which yields information on 
diffusion constants. In contrast with (20| . we consider an 
out of equilibrium situation, because our observable: the 
coordinate of the particle does not reach an equilibrium, 
since the system is infinite. 



II. STOCHASTIC MODELS 



Fractional Brownian Motion 



Fractional Brownian motion is generated from frac- 
tional Gaussian noise, like Brownian motion from white 
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noise. Mandelbrot and van Ness |al defined fBM with 



BH{t) 



[it 



{-T)"-^]dB{T) 



(1) 



where T represents the Gamma function and < H < 
1 is called the Hurst parameter. The integrator B 
is ordinary Brownian motion. Note that B is recov- 
ered when taking H = 1/2. The right hand side of 
Eq. ID) is the sum of two independent Gaussian pro- 
cesses. In the definition, for the first Gaussian pro- 
cess, we identify the so-called fBM of Riemann-Liouville 
type E2f>l. Standard fBM, i.e. Eq. is the only 

Gaussian self-similar process with stationary increments 
The variance of Bnit) is 2DHt^", where Dh ^ 
(r(l - 2H) cos{Htt))/{2Htt). In the following, for some 
given H, we denote the trajectory sample of fBM x{t). 
The properties that uniquely characterize the fBM can be 
summarized as follows: x(t) has stationary increments; 
a;(0) = and {x{t)) = for t > 0; {x^{t)) = 2DHt^" for 
t > 0; x{t) has a Gaussian distribution for t > 0. From 
the above properties, the covariance fmiction is |27i] . 



{xih)xit2))=DHitr +4" -\tl-t2r), tl,t2>0. 

(2) 

The non-independent increment process of fBM, called 
fractional Gaussian noise (fGn), is given by 



dx{t) 
dt ' 



t > 0, 



(3) 



which is a stationary Gaussian process and has a stan- 
dard normal distribution for any t > 0. The mean 
(^(i)) = and the covariance is 



= 2DhH{2H - - t2| 



2H-2 



tlM > 0. 

(4) 



£_{t) is fGn defined in Eqs. |3]) and g]), 1/2 < < 1 is the 
Hurst parameter, 7 > is a generalized friction constant. 
Eq. (5) is called a fractional Langevin equation since 
the memory kernel yields a fractional derivative of the 
velocity (hint use Laplace transform or see P,[2|). Note 
that ii < H < 1/2 the integral over the memory kernel 
diverges, and hence it is assumed that 1/2 < _ff < 1. This 
leads to sub-diffusive behavior (y^) ~ t^^^^ . An over- 
damped fractional Langevin equation, where Newton's 
acceleration term is neglected reads [1| 



0^-lj\t-rr"-'^dr + r^-m- 



(6) 



In what follows we investigate ergodic properties of the 
processes x{t),y{t) and z{t). 



III. ERGODIC PROPERTIES 



We consider the time average 



J^-^[x{t' + A)-x{t'rdt' 



SHx{t)) = 



t 



(7) 



which is a random variable depending on the stochastic 
path x{t). Here A is called the lag time. As is well- 
known, for normal Brownian motion, the ensemble av- 
erage mean square displacement is = 2Dit while 
the time average mean square displacement of the single 
trajectory 5'^{x{t)) = 2Z?iA in statistical sense and in 
the long measurement time limit. Hence we may use in 
experiment a single trajectory of a Brownian particle to 
estimate the diffusion constant Di. Will similar ergodic 
behavior be found also for fBM? 



Fractional Brownian Motion 



B. Fractional Langevin Equation 

The standard Langevin equation with white noise can 
be extended to a generalized Langevin equation with a 
power-law memory kernel. Such an approach was re- 
cently used to model dynamics of single proteins by the 
group of Xic [8] and can be derived from the Kac-Zwanzig 
model of Brownian motion [28| . The under-damped frac- 
tional Langevin equation reads 

where according to the ffuctuation dissipation theorem 




2DhH{2H -ly 



If the average of the random variable (5^ is equal to the 
ensemble average (x^), and if the variance of 5'^ tends to 
zero when the measurement time is long the process is 
ergodic, since then the distribution of 5'^ tends to a delta 
function centered on the ensemble average. For fBM, 
using Eqs. © and ID) 

hence ((5^) = {x"^) for all times. The variance of 5'^{x{t)) 
is, 

Var(5^(x(i))) - - (^{x{t))) . (9) 



3 




FIG. 1: The function of k{H) Eq. (fT2)) . 



A dimensionless measure of ergodicity breaking (EB) 
is the parameter, 



S^x{t)) 



(10) 



which is zero in the hmit t — > oo if the process is ergodic. 
In the next sub-section we derive our main result, vaUd 
for large i, we find: 



EB(x(i)) 



0<H<1 



fc(i?)flnt, H=l 



(11) 





FIG. 2: The EB parameter Eq. (fTTjl for fractional Brownian 
motion x{t) versus time t, for different values of Hurst ex- 
ponent. Here we present exact results obtained directly by 
calculating Eqs. (|8ll5ll6p . The two solid lines (red on line) 
are the asymptotic theory Eq. (|lip for H = 0.2 and H = 0.9. 
In the ballistic limit _H" — s- 1 we get non-ergodic behavior. No- 
tice that for H < 3/4, and long t the curves are parallel to 
each other due to the EB ~ law valid for H < 3/4, while 
for H > 3/4 the slopes are changing as we vary H . The lag 
time is A = 1. 



FIG. 3: We simulate IBM and present the time average 
5'^{x{t)) versus A (dotted curves blue on line). We show 23 
trajectories, the solid line in the middle being the average 
of the trajectories. We observe 5'^{x{t)) oc A^''*, similar to 
that in [ly, |Tl|]. The measurement time is t = 10*, H = 
3/8, Dh ~ 1/2. Line with a slope of 0.75 is drawn to guide 

the eye. We also show {5^) ± \Jvax{P) (two solid lines red 
on line) obtained from Eqs. ((8] I23|l which give an analytical 
estimate on the scatter of the data. 



where 



k{H) 



4i/2(2iJ_i)2^^, 



(12) 



'.4H-3 iH-2> 



1<H<1. 



r 



The most remarkable result is that k{H) diverges when the properties of fractional Brownian motion, see Fig. [T] 
H 3/4, H — 3/4 marks a non smooth transition in Notice that k{H) — > when _ff ^ so the asymptotic 
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convergence is expected to hold only after very long times 
when H is small, since then the diffusion process is very 
slow. When _ff — > 1 we have EB(x(t)) ^ 2 indicating 

ergodicity breaking, see Fig. O 

Fig. [3] displays the simulations of 5'^{x{t)) showing the 
randomness of the time average for finite time measure- 
ments. In this simulation we generate single trajectories 
using the Hosking method 34], and then perform the 
time average to find P. The Fig. mimics the experi- 



mental results on single lipid granule in a yeast cell and 
of mRNA molecules inside a living E-coli cells [13, [lH , 
where H = 3/8 was recorded. Note however that the 
scatter of the experiments data seems larger (see Figs, 
in (lol . [HI), at least with the naked eye. Further we 
did not consider in our simulations the effect of the cell 
boundary. Direct comparison at this stage between ex- 
periments and stochastic theory is impossible, since the 
number of measured trajectories is small. 



I 

B. Derivation of Main Result Ea.(fTT|) 



From Eq. 



2\ f*-^dhf;^-^dt,{[xit,+A)-xit,)nx{h+A)-xit2)?) 



(13) 



Using Eq. ^ and the following formula for Gaussian process with mean zero [sij . 

{x{h)x{t2)x{t3)x{U)) = {x(h)x{t2)){x{t3)x{U)) + {x(U)x{t3)){x{t2)x{U)) + {x{tl)x{U)) {x{t2)x{t3)) , 

we obtain 

{[x{h + A) - x{h)]^[xit2 + A)- xit2)f) = iDj.A^" + 2D%{\ti + A - tap^ + 1^2 + A - h\^" ~ 2\ti - t2\^"y. (14) 
From Eqs. (|8l9ll3ll4p . we have 



2A 



Yav{5\x{t))) = ADjf{ I {t - A ~ t'){{t' + Af" + \t' - A\^" - 2t'^"}^dt' } /{t - A) 



H 



t~A 
2A 



Vi 

{t-A- t'){{t' + A)2^ + {t' - Af" - 2t'^"Ydt' I /{t - Af . 



(15) 



(16) 



V2 



When t >> A we may approximate the upper limit in the integral of V2 with t, and l/(t — A)^ 1/t. We then make 
a change of variables according to x — {t — t') /t and find 



V2=ADlt^" [\{l-x) 
Jo 

We expand in A/t to second order and find 



A 



t{l-x) 



2H 



1 - 



A 



t{l-x) 



2H 



-2 



d.T 



V2 ~ ^Dj,t^" (^j^ ' i/2 (2ij _ 1)2 _ x)^"-^dx. 



(17) 



(18) 



The integral is finite only if > 3/4 hence for ^ 3/4 we will soon use a different approach. We see that V2 ~ ^ 
while it is easy to show that V\ ~ Xjt hence for > 3/4 we find after solving the integral 



Var((52(a;(0)) 16D%t*" { j ] H\2H - 



1 1 



AH -3 AH -2 



Now we write the variance as 



Ya.r{5^{x{t))) 



AD 



2 /.t-A 
H 



(t-A - t') {t' + A) + \t' - A\^" ~ 2t' 



dt'. 



(19) 



(20) 



5 



Changing variables according to r = t'/A we find 



Var(J2(a;(t))) 
The correction term is 



(t - A Jo 



(l + r)^^ + |l-r| 



2H _ 27-2^f 



Corr. 



4r)2 /-VA-i 



{l + rf + ll-Tp" -2t'» 



Taking the upper hmit of the integral in Eq. (pij) to co we find that for iJ < 3/4 and long times 



Yai {6^ {x{t))) r^ADjj A"" [jjj dr 



A 



il + rf" + \l-T\ 



2H 2t^^ 



(21) 



(22) 



(23) 



This is because Corr ~ tmax{4H-4,-2} 

prove this in the following) and this term is smaller than the leading term 

which has a 1/i decay, since H < 3/4. 
Now we estimate the correction term Eq. 



i2 



t/A-l 



dr 



(l + T)'^ + |l-rp« -2t 



2H n9H 



i/A-l 



dr 



2H o^2H 



Using the Lagrange reminder of Taylor expansion in 1/t, when iJ 5C | we have 



t/A-l 



dr 



(l + T)'^ + |l-r|^^ -2t 



2H o^2H 



T. 

(24) 



1 



t/A-l 



I /-t/A-l 



dr 



1 

1 + - 

T 



2ff 



1 

1 - - 

r 



^4ff+l 



^max{4ff-4,-2} 



dr [(1 + 0'^-' + (1 - 0'"-'] - ^ e [0, ^] 



(25) 



For H = J we use Eq. (|2ip . however now we expand to third order and find 



Yai{S^{x{t))) - 16L'|fi?2(2H- l)2A-'lnt - 



while the correction term Corr ~ t ^ [see Eq. (^5])] is negligible. Using Eqs. (|19l23l26p we derive Eq. pT|) . 



(26) 



C. Over Damped Fractional Langevin Equation Then since (^(t)} = we have {z{t)) — 0, and 

We now analyze the over-damped fractional Langevin {z(ti)z{t2)) — Dpiti + ^2 — |^i— ^2! ), (29) 
Eq. ([6]), we can rewrite it in a convenient way as 



where 



where D = d/dt, and D^^^^ is the Riemann-Liouville 
fractional integral of 2H — 1 order. Using the tools of 
fractional calculus [2§|, we get 



Dp = 



ksTi: csc[t:{2 - 2H)] 



7(2 - 2H)V^{2H - l)r2(2 - 2H) ' 



^T{2H - l)z{t) = 77 ■ D^"-^^{t) 



From Eq. ([29]) we learn that Eq. ((27)) exhibits the same 
behavior as fBM Eq. ^ in the sub-diffusion case. Note 
that for fBM (a;^) ^ t° with a = 2H while for the frac- 
(28) tional Langevin equation {z"^) ~ with a = 2 — 2H, 
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lo' 




10"' lo' lo' 10' 

A 

FIG. 4: The time average 5^{y{t)) is a random variable de- 
pending on the underlying trajectory. Total 23 trajecto- 
ries, besides the solid line with mark o being the average of 
the 23 trajectories, are plotted. The measurement time is 
t = 10^,2 ~2H ^ 0.75, Dh = 1/2, fcs = 1, m = 1, 7 = 
1, fo = 1, r = 1. Lines with slopes of 2.0 (ballistic mo- 
tion in short times) and 0.75 (sub-diffusion for long times) 
are drawn to guide the eye. For long A the behavior of the 
under-damped motion is similar to usual fractional Brownian 
motion. 



and of course the diffusion constants have different de- 
pendencies on parameters of the noise. However these 
minor modifications do not change our main result for 
< H < 1/2 obtained in the previous section (only 
switch the value 2H to 2 — 2H). See this note that the EB 
parameter depends on the behavior of correlation func- 
tion Eq. ^ and the latter are identical for the processes 
x{t) and z{t) in the sub-diffusion case, so EB(z) ^ EB(a;). 



D. Under Damped Fractional Langevin Equation 

We now analyze the fractional Langevin equation with 
power-law kernel, namely, Eq. 

= -lj\t - rf"-'%dr + rj ■ m, (30) 

with dy{0)/dt = vq, y{0) ~ 0, where vq is the initial 
velocity. The solution of the stochastic Eq. ([50)1 is 

y{t) = Vo-t-E2HA~l-t^") 

+ ^ /J(t - r) . E,hA-1 ■ it - r?")i{r)dT, 

where 7 = {^T{2H — l))/rn and the generalized Mittag- 
LefSer function is 

f (an + U) 

and Ea.pi-t) ~ {tT{(3 - a))~^ when t +00. 
We have 

fl-2H 

{y{t)) ^vo-t- e^hA-i ■ t'") - ^ r{2-2Hy (^^^ 

and 

(32) 

2fcj3T ,2-2H 

~ 7r(2H-l)r(3-2H) ' '' ' 

where the thermal initial condition: Vq ~ kBT/m is as- 
sumed. 

Note that for short times we have (if'it)) ~ 
{kBT/m)f. Eqs. §^ and were found ^,^]. 
The covariance function of y{t) reads 



{y{ti)y{h)) = viht2E2HA-ltl")E2HA-ltl") 

drds ■ (h - T)E2H.2{-liti ~ rf")it2 - s)E2Ha{-l{t2 ~ s)2«)|r - s\^"-^ 



When ti, t2 tend to infinity. 



^y^'^^y^'^^^ - ,rH2H-T)TH2-2H) i I • - ^^"'"^'^ ~ - ^1'""^' (^4) 



I 

i.e., the covariance of y{t) approximates to the ones of 

z(t), so we can expect in the long-time limit s\ ,79/ -,\ x> 



lo' lo' 10' 

t 

FIG. 5: The ergodicity breaking parameter EB(j/) versus t. 
Simulations of 200 trajectories were used with 2 — 2H — 
0.75, A = 10, fcs = 1, m = 1, 7 = 1, vq = 1, T = 1. The 
stars * is the theoretical result Eq. with corresponding 

parameter values (without fitting). 



and 

EB(y) ~ EB(z) ~ EB(a;). (36) 

The simulations [s^, Fig. [4l confirms Eq. p5|) and Figs. 
[5]and[6l support Eq. ([36|) . Note that for short times we 
have a baUistic behavior for y{t) (see Fig. 4), but not for 
z{t) and x{t), so clearly both A and t must be large for 
Eq. ^ to hold. 



IV. DISCUSSION 

We showed that the fractional processes x{t), y{t) 
and z(t) are ergodic. The ergodicity breaking parame- 
ter decays as a power-law to zero. In the ballistic limit 
(x^) ^ non ergodicity is found. For the opposite local- 
ization limit {x^) ~ t° (i.e. ^ for fBM) the asymp- 
totic convergence is reached only after very long times. 
Our most surprising result is that the transition between 
the localization limit and ballistic limit is not smooth. 
When iJ = 3/4 the EB is changed and the amplitude 
k{H) diverges. Other very different critical exponents of 
fractional Langevin equations were recently found in 
There the critical exponents mark transitions between 
over-damped and under-damped motion. So stochastic 
fractional processes possess a zoo of critical exponents. 
Let us now compare between our results and those de- 
rived based on the CTRW model [l^ • The most striking 



10"' 10' 10' 10' 

A 



FIG. 6: The ergodicity breaking parameter EB(j/) versus A. 
Total 200 trajectories are used to compute the average and 
variance, the measurement time t = 10*, 2 — 2H — 0.75, fcs = 
1, m = 1, 7 = 1, Wo = 1, T — 1. The stars * are the theo- 
retical result Eq. (|lip with corresponding parameter values. 
We see that results found from the under-damped Langevin 
equation converge to our analytical theory based on fBM. 



difference is that for the CTRW model we have non er- 
godic behavior even in the long time limit. It is attempt- 
ing to conclude that this indicates that the underlying 
stochastic motion for the mentioned experiments in the 
cell is of CTRW nature. However, as mentioned in the 
introduction experiments are conducted for finite times, 
and hence what may seem as a deviation from ergodic be- 
havior may actually be a finite time effect. Here we gave 
analytical predictions for the deviations from ergodicity 
for finite time measurement, based on three fractional 
models. The EB parameter depends on measurement 
time and lag time, and can be used to compare experi- 
mental data with predictions of fractional equations (the 
EB parameter for the CTRW is given in [IJ). It should 
be noted however that for sub-diffusion in the cell, effects 
of the boundary of the cell, may be important, and these 
effects where not considered in this text. Another impor- 
tant difference is that for an infinite system we have for 
the CTRW (6^) ~ A/P, so the time average procedure 
yields a linear dependence on A and an aging effect with 
respect to the measurement time. Hence for CTRW an 
anomalous diffusion process may seem normal with re- 
spect to A [13, [13, lla, [131 • In contrast for the fractional 
models we investigated here, we have (5^) ~ A" which is 
the same as the ensemble average (x^) ^ t". The main 
difference between the two approaches, is that the CTRW 
process is non-stationary. It would be interesting to in- 
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